### Make stylized damage figures (Figure 3)
dmg.t.grid = c(1, 25, 50,100,150,200,250,300)
pulse.pattern = exp(cumsum(c(0.05, 0.08, 0.04, 0, 0, 0, 0, 0)*c(dmg.t.grid[1], diff(dmg.t.grid))))
period.length = c(1, diff(dmg.t.grid))
y0 = 0.05 # just scales things up and down. tweak to get a "reasonable" stylized SCC value.

g.examples = c(0, 0.02, 0.04)

stylized.curves = function(beta, rho, eta) {
  md.examples.undisc = matrix(NA, nrow=length(dmg.t.grid), ncol=3)
  md.examples.undisc[,1] = pulse.pattern*y0*exp(beta*g.examples[1]*dmg.t.grid)
  md.examples.undisc[,2] = pulse.pattern*y0*exp(beta*g.examples[2]*dmg.t.grid)
  md.examples.undisc[,3] = pulse.pattern*y0*exp(beta*g.examples[3]*dmg.t.grid)
  
  
  md.examples.disc = matrix(NA, nrow=length(dmg.t.grid), ncol=3)
  md.examples.disc[,1] = exp(-(rho+eta*g.examples[1])*dmg.t.grid)*md.examples.undisc[,1]
  md.examples.disc[,2] = exp(-(rho+eta*g.examples[2])*dmg.t.grid)*md.examples.undisc[,2]
  md.examples.disc[,3] = exp(-(rho+eta*g.examples[3])*dmg.t.grid)*md.examples.undisc[,3]
  
  md.examples.disc.const = matrix(NA, nrow=length(dmg.t.grid), ncol=3)
  md.examples.disc.const[,1] = exp(-0.03*dmg.t.grid)*md.examples.undisc[,1]
  md.examples.disc.const[,2] = exp(-0.03*dmg.t.grid)*md.examples.undisc[,2]
  md.examples.disc.const[,3] = exp(-0.03*dmg.t.grid)*md.examples.undisc[,3]
  
  return(list(md.examples.undisc=md.examples.undisc, 
              md.examples.disc=md.examples.disc, 
              md.examples.disc.const=md.examples.disc.const))
}
beta=1
page.scale = 0.9
pdf(figure.dir%&%'/Discounting_Example_Four_Panel_'%&%today()%&%'.pdf',
    family='serif', width=14*page.scale, height=8.5*page.scale) 

m <- rbind(c(0, 1, 0), c(2, 3, 4))
print(m)
layout(m)

sc.out = stylized.curves(beta=1, rho=0.03 - 1*g.examples[2], eta=1)

par(mai=c(0.6,0.82, 0.25, 0.1)) # BLTR. MArgin in Inches
plot(sc.out$md.examples.undisc[,3] ~ dmg.t.grid, type='n', ylim=c(1e-2, 1e6),
     xlim=c(1,300), cex.axis=1.5, cex.lab=2,
     ylab='Undiscounted future damages ($)',xlab='Years', log='y', yaxt='n', 
     main='')
axis(2, at=10^(-2:6), labels=parse(text='10^'%&%-2:6), tick = T, cex.axis=1.2)
abline(h=10^(-2:6), col='lightgray', lty='dotted')
grid()
lines(sc.out$md.examples.undisc[,3] ~ dmg.t.grid, type='l', col=colorspace::darken(3, 0.5), lwd=2)
lines(sc.out$md.examples.undisc[,2] ~ dmg.t.grid, type='l', col=2, lwd=2)
lines(sc.out$md.examples.undisc[,1] ~ dmg.t.grid, type='l', col=colorspace::lighten(4, 0.5), lwd=2)
legend('topleft', legend=c('4% growth rate','2% growth rate','0% growth rate'), 
       lty=1,   lwd=3, col=c(colorspace::darken(3, 0.5), 2, colorspace::lighten(4, 0.5)), bg='white', cex=1.25)

for (fig.cell in 1:3) {
  if (fig.cell==1) etas.plot = c(0.8, 0.6)
  if (fig.cell==2) etas.plot = 1
  if (fig.cell==3) etas.plot = c(1.2, 1.4)
  
  sc.out = stylized.curves(beta, 0, 0)
  # eta and rho don't matter here since we're just using the function to get the constant discounting cases
  
  par(mai=c(0.6,0.82, 0.25, 0.1)) # BLTR. MArgin in Inches
  plot(sc.out$md.examples.disc.const[,1] ~ dmg.t.grid, type='l', col=colorspace::lighten(4, 0.5), lwd=2, ylim=c(1e-4, 1e2),
       xlim=c(1,300), cex.axis=1.5, cex.lab=2,
       ylab='Discounted future damages ($)',xlab='Years', log='y', yaxt='n')
  axis(2, at=10^(-4:3), labels=parse(text='10^'%&%(-4:3)), tick = T, cex.axis=1.25)
  abline(h=10^(-4:3), col='lightgray', lty='dotted')
  grid()
  lines(sc.out$md.examples.disc.const[,3] ~ dmg.t.grid, type='l', col=colorspace::darken(3, 0.5), lwd=2)
  lines(sc.out$md.examples.disc.const[,2] ~ dmg.t.grid, type='l', col=2, lwd=2)
  lines(sc.out$md.examples.disc.const[,1] ~ dmg.t.grid, type='l', col=colorspace::lighten(4, 0.5), lwd=2)
  
  
  for (eta in etas.plot) {
    line.type = switch(which(eta==etas.plot), 2, 3) # first...dashed, second...dotted
    rho = 0.03 - eta*g.examples[2]
    sc.out = stylized.curves(beta, rho, eta)
    
    lines(sc.out$md.examples.disc[,3] ~ dmg.t.grid, type='l', col=colorspace::darken(3, 0.5), lwd=2, lty=line.type)
    lines(sc.out$md.examples.disc[,2] ~ dmg.t.grid, type='l', col=2, lwd=2, lty=ifelse(eta==1, 0, line.type))
    lines(sc.out$md.examples.disc[,1] ~ dmg.t.grid, type='l', col=colorspace::lighten(4, 0.5), lwd=2, lty=ifelse(eta==1, 4, line.type))
  }
  
  if (all(etas.plot==c(0.8,0.6)))  {
    legend('bottomleft', legend=c('Constant 3% Discounting',
                                  expression(paste('Discounting Rule, ', eta, '=', 0.8)),
                                  expression(paste('Discounting Rule, ', eta, '=', 0.6))),
           lty=1:3, lwd=2, bg='white', cex=1.25)
  }
  if (all(etas.plot==1))  {
    legend('bottomleft', legend=c('Constant 3% Discounting',
                                  expression(paste('Discounting Rule, ', eta, '=', 1))),
           lty=1:2, lwd=2, bg='white', cex=1.25)
  }
  if (all(etas.plot==c(1.2, 1.4)))  {
    legend('bottomleft', legend=c('Constant 3% Discounting',
                                  expression(paste('Discounting Rule, ', eta, '=', 1.2)),
                                  expression(paste('Discounting Rule, ', eta, '=', 1.4))),
           lty=1:3, lwd=2, bg='white', cex=1.25)
  }
}
dev.off()
